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SUMMARY 


The explicit transient dynamics technology in use today for simulating the impact and subsequent 
transient dynamic response of a structure has its origins in the "hydrocodes" dating back to the late 1940's. 
The growth in capability in explicit transient dynamics technology parallels the growth in speed and size of 
digital computers. Computer software for simulating the explicit transient dynamic response of a structure 
is characterized by algorithms that use a large number of small steps. 

In explicit transient dynamics software there is a significant emphasis on speed and simplicity. The 
finite element technology used to generate the spatial discretization of a structure is based on a compromise 
between completeness of the representation for the physical processes modelled and speed in execution. 
That is, since it is expected in every calculation that the deformation will be finite and the material will be 
strained beyond the elastic range, the geometry and the associated gradient operators must be 
reconstructed, as well as complex stress-strain models evaluated at every time step. As a result, finite 
elements derived for explicit transient dynamics software use the simplest and barest constructions 
possible for computational efficiency while retaining an essential representation of the physical behavior. 
The best example of this technology is the four-node bending quadrilateral derived by Belytschko, Lin and 
Tsay (1984). 

Today, the speed, memory capacity and availability of computer hardware allows a number of the 
previously used algorithms to be "improved." That is, it is possible with today's computing hardware to 
modify many of the standard algorithms to improve their representation of the physical process at the 
expense of added complexity and computational effort. 

The purpose of this presentation is to review a number of these algorithms and identify the 
improvements possible. In many instances, both the older, faster version of algorithm and the improved 
and somewhat slower version of the algorithm are found implemented together in software. Specifically, 
the following seven algorithmic items are examined: 

1) The invariant time derivatives of stress used in material models expressed in rate form. 

2) Incremental objectivity and strain used in the numerical integration of the material models. 

3) The use of 1 -point element integration versus mean quadrature. 

4) Shell elements used to represent the behavior of thin structural components. 

5) Beam elements based on stress-resultant plasticity versus cross-section integration. 

6) The fidelity of elastic-plastic material models in their representation of ductile metals. 

7) The use of Courant subcycling to reduce computational effort. 
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INVARIANT TIME DERIVATIVE OF STRESS 


To account for the fact that bodies subjected to large displacements undergo significant rigid body 
rotations, material models in rate form rely on an invariant time derivative of the stress. For example, a 
linear elastic material expressed in terms of rates has the form: 




= C rsmn d 


mn 


( 1 ) 


where 0)to„ is a rotation rate, dmn is the stretching and C rsmn , the material’s elastic modulus. To date the 
majority of simulations have used the Jaumann invariant derivative which uses the spin for the rotation 

rate, COfcm = ( Vjfc.m ~ 

This formulation is very sensitive to shear strains that are greater than 20% in the presence of rotations 
(Dienes, 1979). 

A more accurate invariant time derivative is the Green-Mclnnis derivative which is based on the 

rotation from the polar decomposition of the deformation gradient, co^ =R kf R^ , where F™ = R mn U na 
(Dienes, 1979). 

Computing the uniform shearing of a 1 x 1 coupon demonstrates the difference in predicted behavior 
using kinematic hardening plasticity for these two invariant time derivatives (Fig. 1). 
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(Note: ln(t/t 0 ) = 5 => C = 10) 

Figure 1 - Uniform shearing of a 1 x 1 coupon 



INVARIANT TIME DERIVATIVE OF STRESS (Cont'd.) 


An examination of the shear stress as a function of shear strain produces very interesting results (Fig. 
2). The results based on using the Green-Mclnnis derivative in Fig. 2 are denoted by "Dienes." As can be 
seen, extreme distortions using the Jaumann derivative lead to physically unrealistic stress variations that 
change sign. The monotone increasing curve denoted by "Dienes" exhibits a physically realistic stress- 
strain representation. This anomalous behavior exhibited by the Jaumann derivative must be avoided in a 
simulation if practical results are to be obtained. 

For the sake of reliability, today's software will offer the Green-Mclnnis form of the invariant time 
derivative of stress along side of the Jaumann derivative for material models based on rate formulations. 


UNIFORM SHEARING — KINEMATIC 28 31 30 20-AUG-83 



Figure 2 - Kinematic hardening plasticity predictions of shear 
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INCREMENTAL OBJECTIVITY FOR STRAIN 


In numerical simulations differential expressions are replaced with incremental expressions. 
Incremental expressions for such things as rotation and strain increments can produce errors even for small 
steps if the properties of the differential expression are lost. For example, when the time derivative of an 
orthogonal rotation is approximated, it is very easy to lose the orthogonal properties associated with the 
differential form. 

The same thing can happen when integrals of strain rates are replaced by sums of strain increments, 
particularly when finite strains are involved. Traditionally, a strain increment is obtained from the 

stretching with A = At d mn , d mn = (t) m ,„ + u„, m )/2. 

This approach to generating strain increments can be sensitive to cyclic strain histories in the presence 
of rotations. 

A more accurate, but computationally more costly approach, is to base the strain increment on the 
symmetric stretch tensor U obtained from the polar decomposition of the deformation gradient, 

F™ = R mn U na . This latter approach has been termed strong objectivity by M. Rashid (1992), who 
refers to the traditional approach as weak objectivity. 

Computing the hoop vibration of a rotating ring provides a good example of the contrasting results that 
can come from the use of the faster, classical formulations that are weakly objective and the more accurate 
and costly formulations that are strongly objective. Figure 3 shows an elastic ring given an initial angular 
velocity of 4,000 radians per second. 



Figure 3 - An elastic ring with an initial 4,000 rad/s angular velocity 
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INCREMENTAL OBJECTIVITY FOR STRAIN (Cont’d.) 


The ring rotates approximately 360 degrees in 1.5 milliseconds. During that time thering oscillates 
radially five times (Fig. 4). Of particular interest is the increase in kinetic energy with time (Fig. 4). 
is a closed system to which no additional energy is being added. The increase in kinetic energy reflects the 
accumulation of errors from the strain increments in the classical formulation that is weakly objective in 
spite of the very small steps in the algorithm. 


Rotatir* Ring (Omega - 4,000 rads/sec) HXEL (Daumann, Elastic) FMA-3D V07. 32 31-m 1992 09:50:09 
Strain Based on Emn - (t(n+l)-t(n» x (Vm,n ♦ Vn,m), WeaV Objectivity 



Figure 4 - Kinetic energy response with only weak objectivity 
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INCREMENTAL OBJECTIVITY FOR STRAIN (Cont'd.) 


The same computation using a formulation that is strongly objective in the sense of Rashid shows the 
expected response (Fig. 5); namely, a kinetic energy that oscillates between two limits that are constant in 
time. This latter calculation shows an exemplary energy exchange between kinetic energy and internal 
strain energy as the ring rotates and oscillates radially. 

In spite of the increased expense, today's software should offer as the default a strain increment 
formulation meeting the requirements for strong objectivity. 


Rotating Ring (Onega - 4,000 rada/sec) HXEL (Polar Deconp, Elastic) 1-3UL-1992 08:51:45 
Strain Based on M. Rashid's Definition of Strong Objectivity 



Figure 5 - Kinetic energy response with strong objectivity 
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ELEMENT INTEGRATION 


In the interest of speed, explicit transient dynamic computer programs routinely use constant stress 
elements. In effect, only the uniform strain states of the element are used to represent the structures 
behavior. Strain states and, consequently, stress states that vary over the element are ignored. Such 
approximations result in complex, short wavelength, stress-free deformation patterns that are called 

hourglass modes. , . . , 

The quickest way to obtain a constant stress element is to evaluate the integrals over the area or volume 
of the element defining the gradient and divergence operators with a single integration point; hence, the 
terminology 1 -point integrated elements. Unfortunately, such elements can fail the Iron s patch test. That 
is, a collection of irregularly-shaped elements subjected to a linear motion on the boundary will not 
produce a constant state of strain in the interior and, therefore, the collection of irregularly- shaped elements 
will not produce a constant state of stress. 

A more accurate approach is to evaluate the integrals over the area or volume of the element defining 
the gradient and divergence operators exactly for a constant state of stress, effectively a projection 
calculation. The result is a much more reliable and well-behaved simulation, albeit requiring the execution 
of a greater number of algebra expressions. The mean quadrature elements while still constant strain 

elements, pass the Iron's patch test. . , 

The greatly reduced excitation of the hourglass modes and the greatly improved hourglass control 
offered by the mean quadrature integration over the 1 -point integration virtually renders obsolete the older 
1 -point integration technology. (In the special case of two-dimensional continuum elements, both 

approaches yield the same results.) . . 

Current users of explicit transient dynamic software should only expect to use l -point integrated 
elements in place of elements based on mean quadrature to obtain "quick and dirty results. 
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QUADRILATERAL SHELL ELEMENTS 


In explicit transient dynamics software there is a significant emphasis on speed and simplicity. That is, 
since it is expected in every calculation that the deformation will be finite and the material will be strained 
beyond the elastic range, the geometry and the associated gradient operators must be reconstructed, as well 
as, complex stress-strain models evaluated at every time step. As a result, finite elements derived for 
explicit transient dynamics software use the simplest and barest constructions possible for computational 
efficiency while retaining an essential representation of the physical behavior. The best example of this 
technology is the four-node bending quadrilateral derived by Belytschko, Lin and Tsay (1984). 

The BLT element is based on a constant stress assumption coupled to a flat plate geometric 
approximation of what is usually a warped geometry (the element's geometry is warped when the four 
nodal points do not define a flat surface). In certain situations the BLT element exhibits shortcomings in 
representing the response of a structure. Two examples are the bending of a twisted beam (Example 1), 
and the response of a hemisphere loaded by opposing forces across the hemisphere's diameter (Example 

A four-node bending quadrilateral has been developed that exhibits the expected behavior in these two 
examples. The element has two properties that provide the expected response: 

1) a warping deformation that possesses proper structural stiffness as opposed to being an hourglass 
mode, and 

2) a derivation that is based on the actual geometry of the element as opposed to treating the geometry 
as flat. 

Example 1. The tip-loaded, twisted cantilever beam is an example of a structure with a nonplanar 
geometry that can be nontrivial to reproduce with the finite element technology available today in explicit 
transient dynamic software (Fig. 6). 


B«m with a 90 Degree Spiral Twiel IPlalelG) 



TIKE = 0 OOOOf+OO Magnification » 7 5134 


Figure 6 - Twisted beam with a constant tip load 
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QUADRILATERAL SHELL ELEMENTS (Cont'd.) 


The tip load is applied suddenly at time equal to zero and held constant in magnitude and orientation. 
Based on a length of 12.0 in., a width of 1.1 in., a thickness of 0.32 in. and the selected properties (a 

Young's modulus of 2.9 x 10 7 psi, a Poisson's ratio of 0.22, and a density of 2.5 x 10 4 lbf-sec 2 /in 4 ), 
the static deflection equals 0.005424 in. (MacNeal and Harder, 1985). The fundamental period equals 8.0 
milliseconds. 

As can be seen from Fig. 7, the BLT element with zero hourglass stiffness predicts an unacceptably 
large deflection amplitude and response period; see also the results reported by Belytschko, Wong and 
Chiang (1989). 


Beam with a 90 Degree Spiral Twist (Belytschko-Tsay Shell) 19:52:22 17~9UL-91 



Figure 7 - Beam tip motion prediction using Belytschko-Tsay shell 
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Displ Z at Nodal Point Number 36 (10**-2) 


QUADRILATERAL SHELL ELEMENTS (Cont'd.) 


The technology developed by KEY Associates predicts very nearly the exact amplitude and period 
(Fig. 8). 


Beam with a 90 Degree Spiral Twist (Platei6) 



Figure 8 - Beam tip motion prediction using improved four-node shell 
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QUADRILATERAL SHELL ELEMENTS (Cont’d.) 


Example 2. A hemisphere loaded by two sets of diametrically opposed forces in the plane of the 
equator, separated by 90 degrees and alternating in sign, is a problem in which both bending and 
membrane deformation occur. The loads enter the structure by generating moments including warping or 
twisting moments in the elements adjacent to the loads. 

The loads are applied suddenly at time equal to zero and held constant in magnitude and orientation 
(Fig. 9). 



Figure 9 - Quadrant of a hemisphere loaded across two diameters 
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QUADRILATERAL SHELL ELEMENTS (Cont'd.) 


Based on the hemisphere's radius of 10.0 in., a thickness of 0.04 in. and the selected properties (a 

Young's modulus of 6.825 x 10 7 psi, a Poisson's ratio of 0.3, and a density of 2.5 x 10 -4 lbf-sec 2 /in 4 ), 
the static deflection equals 0.094 inches (MacNeal and Harder, 1985). 

As can be seen from Fig. 10, the BLT element with zero hourglass stiffness predicts an unacceptable 
cyclic accumulation of displacement; see also the results reported by Belytschko, Wong and Chiang 
(1989). 


Diameter Forces On A Hemisphere (BelytschHo-Tsay Shell) 19:56:34 17-JUL-91 



Figure 10 - Radial motion prediction using Belytschko-Tsay shell 
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QUADRILATERAL SHELL ELEMENTS (Cont’d.) 


The technology developed by KEY Associates predicts a steady periodic (constant amplitude) result 
that is within 0.2% (Maximum_Intemal_Energy / 4.0) of the expected deflection (0.094 inches) (Fig. 11). 

Remarks. This element is suitable for both large in-plane and bending strains. These particular 
calculations were carried out without any hourglass control suggesting this element will not be overly 
sensitive to the development of hourglass distortions. In addition, this technology is based on six (6) 
degrees of freedom per nodal point meaning that no further numerical constraints or adaptations are needed 
to eliminate the "drilling" rotation to obtain satisfactory results. Modeling such additional physical features 
as edge beams, or modeling folded plate structures does not present any particular difficulty. 

llie low level of loading in both of these examples results effectively in infinitesimal strains. The 
material response is linear elastic. However, both of these examples require a proper computation of the 
gradient and divergence operators to obtain the correct results. The examples are sensitive indicators of the 
correctness of the representation for the element geometry and the element twisting stiffness. 

While the BLT four-node bending quadrilateral remains a computationally efficient element, there are 
numerous applications for which an accurate representation of the warped geometry and the twisting 
stiffness is essential to obtaining a satisfactory result. Without a doubt, up-to-date software should contain 
an efficient four-node bending quadrilateral with the capabilities discussed here. 


Diameter Forces On A Henisphere <Pl«tel6, 1. 0/2/1) JAWS- 3D V06. 01 16-NOV-1991 06:32:12 



Figure 1 1 - Radial motion prediction (= I E/4) using improved four-node shell 
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BEAM ELEMENTS 


Beam elements, while on the surface appearing to be very elementary structural components, can be 
very intense from a computational standpoint when finite strains and inelastic material behavior are 
modeled. A portion of the computational complexity comes from the need to evaluate the stress-strain 
behavior at many points over the cross section (16 integration stations is a typical number). Membrane, 
bending and torsional stress resultants are produced from weighted sums over these points. 

A very fast and direct method of obtaining stress resultants for elastic-plastic material behavior is to 
introduce stress-resultant plasticity. In effect, the beam cross-section is either completely elastic or 
completely plastic. For simulations in which the beam is deformed well into the plastic range, 
considerable computational effort can be eliminated with stress-resultant plasticity (the evaluation of 
elastic-plastic stress-strain models at numerous individual points over the cross section is not required). 

If only a portion of the beam cross section yields plastically, stress-resultants plasticity will not capture 
any of the detail. 

Quality software seeking to provide options for both rapid results and accurate results will offer both 
beam formulations. 


180 



COURANT SUBCYCLING 


The central difference time integrator is only conditionally stable. That is, the integration procedure 
must be used with a time step that is less than the critical time step. The critical time step is very nearly 
equal to the minimum transient time for a disturbance to cross between any two nodal points in the mesh. 
The algorithm is very effective for an excitation that uses the maximum resolution the mesh can provide; 
for example, a shock wave. 

There are two very common circumstances when an explicit time integration algorithm becomes less 
efficient: first, when there are significant differences in the spacing between nodal points over the mesh 
and, second, when there are significant differences in material stiffness over the mesh. It is not 
uncommon to have differences in critical time steps as much as a hundred to one over a mesh. The 
consequence is that the central difference time integration must function with the smallest critical time step. 

Fortunately, Courant subcycling may be introduced in order to recover much of the efficiency offered 
by explicit algorithms. In Courant subcycling each finite element is integrated with the largest time step 
permitted by the local critical time step. Thus, small, stiff elements are integrated with many small time 
steps, and elsewhere in the mesh large, soft elements are integrated periodically with their larger time 
steps. 

The benefit is significant since the time to perform a calculation drops and the answers are more 
accurate since the central difference algorithm performs as close to the critical step as is practicable 
everywhere. Figure 12 shows an example of a domain in which elements are integrated with two different 
time steps. In this case the ratio is only 2 to 1 , but the principle is demonstrated. 

Clearly, software intended for a wide variety of applications should provide Courant subcycling. 



Figure 12 - Element integration bins for Courant subcycling 
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PLASTICITY IMPLEMENTATIONS 


Traditional elastic-plastic material models are based on a linear representation of plastic hardening and 
radial-return numerical integration. While an exceptionally effective approach for large plastic strains, it is 
somewhat approximate for moderate plastic strains (plastic strains less than a ten times elastic strains. 
Ductile metals can exhibit a very "rounded" stress-strain curve when first yielding plastically. An elastic- 
linear strain hardening model with its transition from elastic to plastic response at a sharp comer does not 
mimic a gradual transition to large plastic strains. 

One should expect a variety of material models in a crash simulation program including the basic linear 
strain hardening model and a model for plastic yielding that has a smooth transition from elastic to plastic 
response. For small plastic strains the models available should include sub-stepping where the time step 
within the material model is broken into smaller steps to control the numerical integration error. 
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SOFTWARE SURVEY 


The following two tables (Tables 1 and 2) are a modest survey of the crash simulation software 
represented by the speakers at this NASA Workshop on Crashworthiness Technology. In brief, the 
software surveyed is as follows: 

1. FMA-3D is a full-feature program available commercially from KEY Associates, 1851 Tramway 
Terrace Loop NE, Albuquerque, NM 87122 that include subcycling. 

2. DYNA-3D is a full-featured program available through a technology sharing agreement from the 
Lawrence Livermore National Laboratory, Attn: Robert Whirley, P.O. Box 808, Livermore, CA 94550. 

3. PRONTO-3D is a limited-feature program available with a license agreement from Sandia National 
Laboratories, Attn: Steven Attaway, Organization 1425, Albuquerque, NM 87185. 

4. SUPER WHAMS is a full-featured program soon to be available commercially from KBS2, Suite 
112, 455 South Frontage Road, Burr Ridge, IL 60521. 

5. DYCAST is a limited-feature program based on an implicit integration algorithm available 
commercially from Grumman Aerospace Corporation, Attn: Allan Pifko, MS A08-35, Bethpage, NY 
11714. 


Table 1. Technology Availability Table 2. Technology Availability 



FMA-3D 

DYNA 

PRONTO 


SUPER WHAMS 

DYCAST 

Time Derivative 




Time Derivative 



Jaumann 

Yes 

Yes 

No 

Jaumann 

Yes 

(no solid 

Green-Mclnnis 

Yes 

No 

Yes 

Green-Mclnnis 

No 

elements) 

Strain Increment 




Strain Increment 



Weak Objectivity 

Yes 

Yes 

No 

Weak Objectivity 

Yes 

(no solid 

Strong Objectivity 

Yes 

No 

Yes 

Strong Objectivity 

No 

elements) 

Element Integration 




Element Integration 



1-Point 

No 

Yes 

No 

1 - Point 

No 

(no solid 

Mean Quadrature 

Yes 

Yes 

Yes 

Mean Quadrature 

Yes 

elements) 

4-Node Shells 




4-Node Shells 



BLT Quad 

Yes 

Yes 

Yes 

BLT Quad 

Yes 

(A l/d&r) 

Improved Quad 

KEY 

EWQ 

No 

Improved Quad 
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Beam Elements 




Beam Elements 



Illyshin Plasticity 

(i/i) 

Yes 

No 

Illyshin Plasticity 

No 

No 

Integrated 

Yes 

Yes 

No 

Integrated 

Yes 

Yes 

Plasticity 




Plasticity 



Linear Hardening 

Yes 

Yes 

Yes 

Linear Hardening 

Yes 

No 

w / Radial Return 




w/ Radial Return 



Smooth Hardening 

Yes 

"Yes" 

No 

Smooth Hardening 

No 

Yes 

w/ Sub-Stepping 




w/ Sub-Stepping 



Time Subcycling 

Yes 

No 

No 

Time Subcycling 

Yes 

Implicit 
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